 /* program to generate theory plots for the bias: Phm/Pmm and Phh/Pmm
  */
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "bias.h"

int main(int argc, char *argv[]) {
  printf("USAGE: bias q ftNL name Np z\n");
  h = 0.7;
  c=299792.458; // in km/s
  k0=0.002/h; // in (h*Mpc^-1)
  kmin=2*M_PI/2400; // k_min_sim in h^-1 Mpc
  /* specify the cosmology */
  Omm=0.27; OmL=0.73; Omb=0.04; sigma8=0.8;
  rhobar=2.775*h*h*1e11*Omm; // M_sun/Mpc^3
  //printf("bias: %g\n", bias(0.05, 0.1, 50000, 8E-10, 3, 1.6, 1));
  bgmin=5.0; dcmin=3.0; chi2min=2000.; chi2t=0.;
  q = atof(argv[1]); ftNL=atof(argv[2]);
  double A=8.20092e-10; double pcmin=0;
  double qfnl;
  qfnl = atof(argv[1])*atof(argv[2]);
  A = calA(qfnl);
   // scan over dc and bg
  printf("%s\n", argv[3]);
  double p;
  p=1.0;
  // put in the value of gaussian bias in p from the measurements already performed. for example to start with work on the ccc, z=1, Np=100 case
  
  for (bg=0.5; bg<15.4; bg=bg+0.01) {
    for (dc = 1.46; dc<1.47; dc=dc+0.01) {
        chi2t =biaschisq(A, dc, q,ftNL, argv[3], bg, p, atoi(argv[4]), atoi(argv[5]));
        if (chi2t<chi2min) {
	        chi2min=chi2t;
	        bgmin=bg;
	        dcmin=dc;
	        pcmin=p;
        }
//	printf("bgmin: %g, dcmin: %g, chi2: %g\n", bg, dc, chi2t); 
    }
  }
  dcmin=1.46;
  dc = dcmin;
  bg = bgmin;
  double bgerror;
  // find the error in bgmin by simply calculating what generates chi^2+1, assume symmetric +-
  chi2t=chi2min;
  bg=bgmin;
  while (chi2t<chi2min+1.0) {
    chi2t = biaschisq(A, dc, q, ftNL, argv[3], bg, p, atoi(argv[4]), atoi(argv[5]));
    bgerror = bg-bgmin;
    bg=bg+0.01;
  }
  printf("dc: %g, bg: %g, bgerror: %g, chi2: %g\n",dcmin, bgmin ,bgerror, chi2min);

  // after getting the best fit generate theory curves and put it in the appropriate directory:
  char fname[150], string[500];
  FILE *fp;
  double k;
  char Opbase[200];
  sprintf(Opbase, "/home/sza5154/Research/feedersim/Theory/Output1024");
  sprintf(fname, "%s/%s/bias/theory.z%d.%d.ratiocm", Opbase, argv[3], atoi(argv[5]), atoi(argv[4]));  
  fp = fopen(fname, "w");
  for (k=kmin; k< 0.08; k=k+kmin/10) {
    fprintf(fp, "%g %g\n", k, bias(k, q, ftNL, A, bgmin, dcmin, p, Mp(atoi(argv[4])),atoi(argv[5])));
  }
  fclose(fp);
  
  
  // generate theory curves for hhbias
  //bgmin=bgmin*1.2; //1.35 is a fudge factor based on inspection of gaussian bias--do a proper chi2 fit to get the factor later.
  sprintf(fname, "%s/%s/bias/theory.z%d.%d.ratiogm", Opbase, argv[3], atoi(argv[5]), atoi(argv[4]));  
  fp = fopen(fname, "w");
  for (k=kmin; k< 0.08; k=k+kmin/10) {
    fprintf(fp, "%g %g\n", k, hhbias(k, q, ftNL, A, bgmin, dcmin, p, Mp(atoi(argv[4])), atoi(argv[5])));
  }
  fclose(fp);
  // generate theory curves for stoch
  // roughly determine the factor:

// search chi2 min varying bg and mu
  dc=dcmin;p=1.0;
  double mu;
  chi2min=1000.;
  for (bg=bgmin; bg<bgmin+1.0; bg=bg+10) {
    for (mu=0.2; mu<1.5; mu=mu+0.01) {
        chi2t =stochchisq(A, dc, q,ftNL, argv[3], bg, p, mu, atoi(argv[4]), atoi(argv[5]));
        if (chi2t<chi2min) {
	        chi2min=chi2t;
	        bgmin=bg;
	        dcmin=dc;
	        pcmin=mu;
	        //printf("dc: %g, bg: %g, mu: %g, chi2: %g\n",dc, bg,mu, chi2t);
      }
    }
  }
  double muerror;
  // find the error in bgmin by simply calculating what generates chi^2+1, assume symmetric +-
  chi2t=chi2min;
  mu=pcmin;
  while (chi2t<chi2min+1.0) {
    chi2t = stochchisq(A, dc, q,ftNL, argv[3], bgmin, p, mu, atoi(argv[4]), atoi(argv[5]));
  //  printf("chi2t: %g, chi2min: %g, mu: %g\n", chi2t, chi2min, mu);
    muerror = mu-pcmin;
    mu=mu+0.01;
  }

  printf("dc: %g, bg: %g, mu: %g, muerror: %g, chi2: %g, NDF: %d\n",dcmin, bgmin, pcmin,muerror,  chi2min, NDFglobal);
  
  sprintf(fname, "%s/%s/bias/theory.z%d.%d.ngstoch", Opbase, argv[3], atoi(argv[5]), atoi(argv[4]));  
  fp = fopen(fname, "w");

  for (k=kmin; k< 0.08; k=k+kmin/10) {
    //printf("k: %g, q: %g, ftNL: %g, stoch: %g\n", k, q, ftNL, stoch(k, q, ftNL, A, bgmin, dcmin, p));
    fprintf(fp, "%g %g\n", k, pcmin*stoch(k, q, ftNL, A, bgmin, dcmin, p, Mp(atoi(argv[4])), atoi(argv[5])));
  }
  fclose(fp);
}

double alpha(double k) { //use k in Mpc^-1 //only case of such; also this is for z=0, the growth(0)=0.76001 factor is included below!
  TFset_parameters(Omm*h*h, Omb/Omm, 2.728);
  //printf("Tranfer(0.001): %g\n", TFfit_onek(k*h, 0, 0));
  return 2*TFfit_onek(k, 0, 0)*k*k*0.76001*c*c/(3*h*h*100*100*Omm); // * hubble means the input k is in h*Mpc^-1!
}


double calA(double qftNL) {
  return 8.20092e-10*1.0/(0.5 + 6.325e-9*sqrt(6.25e15 + 3.2e8*pow(qftNL,2.0)));  // this factor has been checked to agree with the ratio measured directly from IC code;
}


double stochchisq(double A, double dc, double q, double ftNL, char* name, double bg, double p, double mu, int M, int z) {
// mu is the extra multiplicative factor
 char fname[150], string[500], fnamegaus[150];
  FILE *fp;
  FILE *fpgaus;
  double k, bk;
  int ndf;
  char Opbase[200];
  sprintf(Opbase, "/home/sza5154/Research/feedersim/Theory/Output1024");
  sprintf(fname, "%s/%s/bias/power.z%d.%d.ngstoch.avg", Opbase, name, z, M);
  if (name=="gaus" || name=="gaus2") {
     sprintf(fname, "%s/%s/bias/power.z%d.%d.stoch.avg", Opbase, name, z, M);
  }
  fp = fopen(fname, "r");

  float col1, col2, col4, chi2=0.;
  ndf=0;
  while (!feof(fp)) {
    if (fgets(string, 499, fp)) {
      sscanf(string, "%g %g %*g %g", &col1, &col2, &col4);
      if (col1>0.005 && col1< 0.04 && col2>0) {
	k=col1; // in h*Mpc^-1
	bk=mu*(stoch(k, q, ftNL, A, bg, dc,p, Mp(M), z));
	chi2+=pow((bk-col2)/col4,2);
	ndf++;
      }
    }
  }
  fclose(fp);
 // printf("chi2: %g, ndf: %d, chi2/ndf: %g\n", chi2, ndf, chi2/ndf);
  NDFglobal=ndf-1;
  return chi2/(ndf-1.0);
}

double biaschisq(double A, double dc, double q, double ftNL, char* name, double bg, double p, int M, int z) {
  char fname[150], string[500];
  FILE *fp;
  double k, bk;
  int ndf;
  char Opbase[200];
  sprintf(Opbase, "/home/sza5154/Research/feedersim/Theory/Output1024");
  sprintf(fname, "%s/%s/bias/power.z%d.%d.ratiocm.avg", Opbase, name, z, M);
  fp = fopen(fname, "r");
  float col1, col2, col4, chi2=0.;
  ndf=0;
  
  while (!feof(fp)) {
    if (fgets(string, 499, fp)) {
      sscanf(string, "%g %g %*g %g", &col1, &col2, &col4);
      if (col1>0.003 && col1< 0.04 && col2>0) {
	k=col1; // in h*Mpc^-1
	bk=bias(k, q, ftNL, A, bg, dc,p,Mp(M),z);
	chi2+=pow((bk-col2)/col4,2);
	ndf++;
      }
    }
  }
  fclose(fp);
 // printf("chi2: %g, ndf: %d, chi2/ndf: %g\n", chi2, ndf, chi2/ndf);
 NDFglobal=ndf-1;
  return chi2/(ndf-1.0);
}

double I1(double k) {
 if (k>0.0055) {
    return (4.-4.*pow(k/kmin, 0.04) )/(-0.04); // works great for k~>2*k_min
  } else {
    k = k*h;
    return 0.000011652015560618931/k/k - 
 0.006689497124679371/k - 2.660055915097187 + 
 1970.0848555601901*k - 101889*k*k;
  }
}

double growth(double z) {
  if (z==0.) {return 0.76001;} else if (z==0.5) {return 0.59455;} else if (z==1.) {return 0.473345;} else if (z==2.) {return 0.327534;} else {printf("no growth function value!\n"); return 0; }
}

double bias(double k, double q, double ftNL, double A, double bg, double dc, double p, double M, int z) {
  // return the scale dependent bias given the parameters
  double qf;
  double facN;
  double ak;
  facN=2.0*A*I2R2(M);
  qf = q/(1.0+A*pow(q*ftNL,2.0)*facN); // check the numerical factor
  ak=alpha(k*h)*growth(z*1.0)/growth(0.);
  return bg+2.0*dc*(bg-p)*qf*q*ftNL*(1+q*ftNL*ftNL*Pvar(k,A)*I1(k))/(1+q*q*ftNL*ftNL*I1(k)*Pvar(k,A))/ak;
  //return bg+2.0*dc*(bg-p)*qf*q*ftNL/alpha(k);
}

double hhbias(double k, double q, double ftNL, double A, double bg, double dc, double p, double M, int z) {
  // return the hh bias i.e. Phh(k)/Pmm(k), then stoch=(bias)^2/hhbias
  double qf, beta, ak;

  ak=alpha(k*h)*growth(z*1.0)/growth(0.);
  double facN;
  facN=2.0*A*I2R2(M);
  qf=q/(1.0+A*pow(q*ftNL,2.0)*facN); // 15.1 is for Mass[8] which works well for halos with particles 100-200; not for all.
  beta=dc*(bg-p);
  return bg*bg+4.0*beta*qf*q*ftNL*(bg/ak+beta*qf*ftNL/ak/ak)*((1.0+q*ftNL*ftNL*Pvar(k,A)*I1(k))/(1.0+q*q*ftNL*ftNL*Pvar(k,A)*I1(k)));
  //return bg*bg + 2.0*beta*qf*q*ftNL*(bg/ak+2*beta*qf*ftNL/ak/ak);
}

double stoch(double k, double q, double ftNL, double A, double bg, double dc, double p, double M, int z) {
  return pow(bias(k,q,ftNL,A,bg,dc,p, M,z),2.0)/hhbias(k,q,ftNL,A,bg,dc,p, M, z);
}

double Mp(int M) {
if (M==100) {return 150.*9.648e11;} else if (M==200) {return 300.*9.648e11;} else if (M==50) {return 75.*9.648e11;} else {return 142.*9.648e11;}
}

double Pvar(double k, double A) {
  return A*pow(k/k0,-0.04);
}

double I2R2(double M) {
  return pow(M, -3.13172)*exp(271.811+64250.2/pow(log(M),2)-6746.78/log(M));
}
